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Abstract 

in  the  afterbody  region  of  a  supersonic  vehicle 
with  propulsive  jets  is  difficult  to  predict  accurately  partly 
because  the  complex  pattern  of  shocks,  expansion  fans, 
and  free  shear  layers  is  difficult  to  resolve,  and  partly 
because  the  turbulence  is  difficult  to  model.  Th^ 
.^'^obf^Vof  numerical  resolution,  and  turbulence  modeling 
are  addressed,  and  'U^isPshowK  that,  while  insufficient 
resolution  causes  some  errors,  larger  errors  can  be 
attributed  to  the  turbulence  model.  The  standard  k-C. 
model,  gives  good  results  for  one  set  of  experimental 
data/despite  its  inability  to  model  many  effects  on  the 
turbulence  which  may  be  present  in  the  base  region.  A 
modification  for  Mach  number  effects  is  derived  and 
verified  against  experimental  data  for  plane  shear  layers. 
When  the  modified  model  is  used  Wv-  the  computation  of 
^-''',the  base  flow,  the  predicted  base  pressure  rises 
significantly.  This  suggests  that  good  results  with  the 
standard  k-^'  model,^  or  any  model  which  does  not  predict 
Mach  number  effectsT^are  due  to  fortuitous  error 
cancellation. 


Introduction 


A  simple  bluff  base  flow  consists  of  a  reversed  flow 
region  adjacent  to  the  base  which  is  bounded  by  a 
propulsive  jet  and  the  freestream,  as  sketched  in  Fig.  1. 
The  base  pressure  and  the  shape  of  the  reversed  flow 
region  are  determined  by  several  competing  mechanisms. 
Entrainment  into  the  freestream  and  propulsive  jet  pumps 
fluid  out  of  the  reversed  flow  region.  The  outer  shear 
layer,  which  emanates  from  the  “shoulder,"  is 
moderately  thick  and  turbulent.  The  inner  shear  layer, 
which  originates  at  the  nozzle  lip,  is  initially  very  thin 
and  probably  laminar,  at  least  in  non-combusting 
laboratory  experiments.  Often  there  is  greater  mixing  in 
the  outer  shear  layer  which  induces  a  preferred  rotation 
in  the  base  region  (clockwise  in  the  sense  of  Fig.  1). 
The  flow  in  the  base  region  is  at  relatively  low  speed, 
and  thus  this  region  is  roughly  at  uniform  pressure. 
Expansion  fans  emanating  from  the  shoulder  of  the  base 
and  the  lip  of  the  nozzle  cause  the  base  pressure  to  be 
substantially  below  the  freestream,  for  the  configuration 
shown.  This  produces  base  drag,  which  can  be  an 
appreciable  part  of  a  vehicle's  total  drag.  The  low 
pressure  adjacent  to  the  base  pulls  fluid  from  the  region 
where  the  shear  layers  merge  and  the  pressure  is  slightly 
higher,  to  make  up  for  fluid  lost  to  entrainment.  Base 
flows  result  from  a  rather  complex  balance  between  these 
competing  mechanisms,  all  of  which  must  be  computed 
accurately  if  good  predictions  are  to  be  obtained. 

The  ability  of  contemporary  Navier-Stokes  methods  to 
predict  this  type  of  flow  was  investigated  by  Petrie  and 
Walker.’  who  ^mpared  true  predictive  calculations  by 
several  workers  '  to  experimental  data.  Not  only  did 
the  predictions  give  poor  agreement  with  the  experiment. 
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but  also  there  were  significant  variations  in  the  results 
from  nominally  similar  prediction  methods.  All  methods 
which  used  the  Baldv.'in-Lomax^  model  underpredicted  the 
base  pressure,  and  the  one  method  which  used  the  k-c 
or  k-W  models  overpredicted  the  base  pressure.  The 
importance  of  good  grid  ^^stering  was  demonstrated. 
More  recent  computations^*  ’  of  the  low-pressure-ratio 
experiment  of  Petrie  and  Walker  have  focused  on  grid 
resolution,  while  still  using  the  Baldwin-Lomax  model, 
and  have  given  improved  results  for  the  base  pressure. 
However,  significant  grid  dependence  was  seen  in  the 
solutions,  and  thus,  the  predictive  capability  of  the 
mathematical  model  cannot  be  evaluated.  At  this  time 
the  ability  to  predict  accurately  the  type  of  base  flow 
sketched  in  Fig.  1  has  not  been  demonstrated. 

Putnam  and  Bissinger  reviewed  predictions  for  a  wide 
range  of  nozzle  afterbody  flows.  They  observed  that 
Navier-Stokes  methods  gave  poor  results  downstream  of 
separation.  Thus,  bluff  base  flows  and  the  associated 
drag  were  not  predicted  very  well.  One  of  their 
conclusions  was  that  Navier-Stokes  methods  are  not 
reliable  for  the  prediction  of  afterbody  drag. 

The  difficulties  in  obtaining  accurate  numerical 
predictions  fall  into  two  general  categories.  The  first  is 
the  problem  of  computing  accurate  solutions  to  the 
partial  differential  equations  (PDEs)  used  to  model  the 
physics,  which,  in  the  present  work,  are  the  compressible 
Reynolds-averaged  Navier-Stokes  equations  with  a  k-e 
turbulence  model.  The  second  area  concerns  the 
accuracy  with  which  the  mathematical  representation 
models  the  physics.  Particular  interest  focuses  on  the 
turbulence  model.  Within  these  two  general  areas,  there 
are  several  aspects  relevant  to  base  flows. 

The  problem  of  obtaining  accurate  solutions  to  the 
model  equations  centers  on  several  narrow  regions  (at 
shocks,  expansion  corners,  and  in  free  shear  layers)  where 
the  solution  varies  rapidly  in  an  otherwise  smooth  flow. 
Strong  normal  and  oblique  shocks  must  be  computed 
without  excessive  numerical  smearing  or  oscillations.  Pre¬ 
shock  Mach  numbers  in  excess  of  five  are  seen  in  the 
present  work,  and  even  higher  Mach  numbers  will  be 
present  in  other  applications.  The  thin  shear  layers 
bordering  the  recirculating  flow  region  need  to  be  resolved 
so  that  the  mean  flow  and  the  turbulence  quantities 
(when  a  multi-equation  model  is  used)  are  unaffected  by 
numerical  diffusion.  The  expansions  at  the  sharp  corners 
of  the  base  often  cause  numerical  difficulties  and  may 
give  rise  to  unusual  and  unexpected  separation 
characteristics.  Turbulence  in  the  base  region  is  difficult 
to  measure,  and  it  is  poorly  understood.  It  is  affected 
by  pressure  gradients,  including  interactions  with  shocks 
and  expansions,  streamline  curvature,  the  merging  of 
shear  layers,  and  high  Mach  numbers.  Individually,  these 
effects  are  somewhat  understood  yet  still  difficult  to 
model:  together  they  may  produce  entirely  unexpected 
behavior.  Present  models  cannot  be  relied  on  to  give 
reasonably  accurate  predictions  over  a  wide  ranges  of 
conditions. 


■  ’.Jl 


Progress  has  been  made  in  some  of  these  areas. 
Differencing  schemes  with  good  shock  capturing  abilities 
have  been  developed  (e.g.  Ref.  10,11).  Adaptive  grid 
methods  have  demonstrated  the  ability  to  give  good 
clustering  in  r^ions  of  high  gradients  (e.g.  Ref.  2).  The 
area  which  will  likely  be  the  most  difficult  to  solve  to 
the  point  at  which  reliable  accurate  engineering 
calculations  can  be  made  is  that  of  turbulence  modeling. 

It  should  be  noted  that  accurate  measurements  in  the 
base  region  are  also  difficult  to  obtain.  Because  of  the 
high  Mach  numbers,  there  are  density  and  temperature 
fluctuations  which  reduce  the  accuracy  of  hot  wire 
anemometry.  Laser  based  instruments  may  also  be 
affected  by  fluctuations  in  the  optical  density,  as  well  as 
particle  seeding  and  particle  lag  problems,  especially  in 
the  high  speed,  low  density  regions  of  the  plume. 
Shadowgraph-type  visualizations,  surface  oil  flows,  and 
surface  pressures  probably  yield  the  most  reliable  data 
from  these  flows  at  present,  despite  their  limitations. 

The  objectives  of  the  present  work  are  to  evaluate 
the  dependence  rf  solutions  on  the  numerical  grid  used, 
to  determine  what  must  be  done  to  eliminate  any 
significant  grid  dependence,  and  to  assess  and  improve 
upon  the  accuracy  of  turbulence  modeling  for  the  base 
region. 

In  the  remainder  of  this  report  the  following  will  be 
described:  the  numerical  method,  the  turbulence  model, 
including  a  modification  for  Mach  number  effects,  and  the 
grid  generation  method.  Then  the  results  of  several 
calculations  will  be  discussed. 


Navier-Stokes  Algorithm 

A  finite  volume  algorithm,  with  a  characteristic-based 
method  of  computing  fluxes  and  an  implicit  time  advance 
method,  is  used  to  solve  the  compressible  Reynolds- 
averaged  Navier-Stokes  equations  and  the  k-e  turbulence 
model  equations.  A  brief  description  of  the  distinguishing 
aspects  of  the  method  is  given  below.  The  algorithm 
and  code  used  in  the  present  work  were  developed  by 
Coakley,  and  the  reader  is  referred  to  Ref.  12  for  details 
of  the  method. 

Base  flows  often  contain  complex  patterns  of  shocks, 
and  it  is  desirable  to  resolve  these  sharply,  without 
excessive  numerical  diffusion^  ^^There  is  a  broad  class  of 
algorithms  for  inviscid  flow,^  ‘  '  often  referred  to  as  high 
resolution  methods,  which  have  very  good  shock 
capturing  abilities.  Strong  shocks  can  be  captured  over 
two  or  three  grid  points  with  little  or  no  oscillation  in 
the  solution,  which  removes  the  need  to  do  case  specific 
grid  clustering  to  obtain  well  defined  captured  shocks. 
Coakley's  method  for  inviscid  flux  calculations  falls  into 
this  class,  and  it  uses  upwind  characteristic  differencing, 
with  first,  second,  or  third  order  spatial  accuracy. 
Optionally,  second  order  central  differencing  can  be  used. 
Except  where  noted,  the  second  order  upwind  method  is 
used  for  inviscid  fluxes.  Viscous  terms  are  always 
computed  with  second  order  central  differencing. 

The  second  order  upwind  difference  scheme  is 
described  here.  For  this  discussion,  the  equations  are 
linearized  so  that  a  brief  yet  clear  explanation  of  the 
differencing  scheme  can  be  given.  The  one-dimensional 
Euler  equations  can  be  expressed  as 

Qt  ♦  ‘  °  (1) 

F  ■=  [/»u,/»u^+p,(/>e+p)u] 


in  which  p  is  density,^  u  is  velocity,  e  is  specific  ir)ternal 
energy,  p  =  (7-l)(e-u72)  is  pressure,  and  7  is  the  ratio 
of  specific  heats. 

The  concern  is  with  the  approximation  of  F*.  The 
flux  vector  F  can  be  split  into  components,  and  F  . 
associated  with  signals  moving  in  the  +x  and  -x 
directions,  with  the  property  F  =  F"*"  +  F  .  The  split 
is  accomplished  by  expressing  F  as  AQ,  in  which 
A  =  8F/dQ,  and  by  performing  eigenvector  manipulations 
on  A.  (The  linear  assumption  means  that  A  is  spatially 
invariant.)  F*  is,  in  effect,  computed  by  conventional 
second  order  upwind  differences  applied  to  F'*'  and  F  . 

In  actual  application  in  the  finite  volume  method  for 
the  nonlinear  case,  the  differencing  scheme  has  a  very 
different  appearance,  and  the  details  of  the  application 
are  crucial  to  the  methods  properties.  Steady  solutions 
of  the  Navier-Stokes  equations  are  obtained  by  time 
asymptotic  integration  with  a  factored  implicit  method. 
A  spatially  variable  time  step  is  used.  Again,  the  reader 
is  referred  to  Ref.  12  for  details  on  the  method. 

Turbulence  Model 

One  possible  approach  to  modeling  the  turbulence  in 
the  base  region  is  to  adopt  the  most  "advanced"  model 
which  is  practical  to  use  and  hope  that  it  captures  the 
critical  turbulence  mechanisms.  The  Reynolds-stress 
transport  model  (RSTM)  is  an  example  of  an  advanced 
model  which  will  capture  some  complex  turbulence 
mechanisms,  and  for  certain  types  of  flows,  a  RSTM  will 
likely  give  very  good  results.  However,  it  can  be  inferred 
from  Ref.  13  that  the  ’’standard’’  RSTM  does  not 
predict  high  Mach  number  effects  nor  does  it  predict 
accurately  the  curving  and  merging  shear  layers,  for 
example.  Ad  hoc  modifications  would  be  required  to 
predict  these  and  other  phenomena  which  are  expected. 
Ad  hoc  modifications  can  also  be  used  to  improve  a  less 
complex  model,  which  may  then  be  more  practical  and 
less  expensive  to  use  than  the  advanced  model.  Either 
approach  introduces  additional  adjustable  coefficients,  and 
neither  approach  has  shown  clearly  superior  results  for 
the  types  of  flows  which  are  of  interest  here.  Thus,  the 
approach  taken  here  is  to  use  the  k-c  model  as  the 
basis  upon  which  to  develop  an  improved  model. 

Two  turbulence  models  are  used  in  the  present  work: 
the  "standard”  high  Reynolds  number  k-E  moder*  and 
the  same  model  with  a  modification  to  account  for  the 
effects  of  Mach  number.  The  former  is  not  presented 
here:  only  the  modification  for  high  Mach  numbers  is 
covered.  There  are  effects  other  than  those  due  to  Mach 
number  which  should  also  be  included  in  a  model  for 
base  flows.  These  will  be  studied  in  the  future. 

The  starting  point  for  the  high  Mach  number 
modification  is  the  RSTM  of  Bonnet  (Ref.  13,  pp.  1408- 
1410).  Bonnet  derived  an  expression  for  the 
compressible,  return-to-isotropy  part  of  the  pressure-strain 
term,  which  is  similar  to  its  incompressible  analogue.  By 
assuming  pressure  fluctuations  are  small  relative  to 
density  fluctuations,  the  model  for  the  complete  pressure- 
strain  term  (compressible  and  incompressible  parts)  is 
reduced  to  the  incompressible  form  with  a  scale  factor 
which  depends  on  the  mean  Mach  number.  On  the 
strength  of  plausible  arguments,  the  scale  factor  is  made 
to  depend  on  the  Mach  number  of  the  fluctuating 
velocity,  u'u'/a,  where  a  is  the  acoustic  speed. 


However,  the ,  algebraic  form  of  the  scale  factor  is 
unchaitged.  The  scale  factor  has  one  coefficient,  which 
was  optimized  by  numerical  experiment,  and  then  the 
model  gave  good  results  for  the  spreading  rate  of  a 
planar  free  shear  layer  up  to  M  =  7.  For  M  >  7 
Bonnet's  numerical  method  was  not  stable. 

In  the  present  work.  Bonnet’s  model  is  reduced  to  a 
form  compatible  with  an  eddy-viscosity  model,  such  as 
the  k-c  model.  To  do  this,  it  is  assumed  that  the  flow 

is _ subject  to  homogeneous  shea'  and  that 

u'u‘/k  =  0.96,  which  is  an  average  of  this  parameter 
for  uniform  shear  flow  and  isotropic  turbulence  in 
incompressible  flow.  With  these  assumptions  and  the 
RSTM  coefficients  Cj  =  1.8  and  C2  =  0.6  (see  e.g.  Ref. 
13)  in  Bonnet's  model,  the  compressibility  modification 
used  in  the  present  work  is  determined.  It  is  a 
modification  of  C^, 

(1  +  C^iW/a2)(l  ♦  (3) 

in  which  Cd  =  8.4  and  C(2  =  6.5.  As  with  with  the 
standard  k-c  i^odel,  Cuq  =  0.09,  and  the  eddy-viscosity 
is  i/i  =  C^pk^/e.  The  variation  of  C^/Cuq  with  Vk/a 
is  given  in  Fig.  2.  The  accuracy  of  Eq.  (3)  is  evaluated 
through  calculations  of  the  spreading  rate  of  a  plane  free 
shear  layer,  which  are  discussed  in  the  Results  section. 
It  should  be  emphasized  that  the  only  numerical 
optimization  done  in  the  development  of  Eq.  (3)  was  that 
done  by  Bonnet.  It  is  suspected  that  Eq.  (3)  may 
underpredict  shear  stresses  at  very  high  Mach  numbers, 
but  this  was  not  studied  in  the  present  work.  Future 
work  may  include  optimizing  and  Cc2< 

Also  shown  in  Fig.  2  is  the  compressibility  correction 
used  by  Pergament  et  al.,‘  which  was  derived  on  a 
purely  empirical  basis.  Their  compressibility  correction 
differs  from  the  present  one  in  two  major  ways.  It 
predicts  a  more  rapid  onset  of  Mach  number  effects. 
This  is  at  least  partially  due  to  the  assumption  in  the 
present  work  about  the  ratio  of  u'u'/k.  If  u'u'/k 
had  been  assumed  to  be  that  for  uniform  shear  flow, 
then  C^/C^O  would  decrease  more  rapidly  with  increasing 
Mach  number.  Secondly,  their  method  bases  the 
correction  on  k,„„f  the  maximum  value  of  turbulence 
kinetic  energy  at  an  axial  station,  whereas  the  present 
method  uses  the  local  value  of  turbulence  energy. 

Boundary  Conditions 

Boundary  conditions  must  be  supplied  at  inflow 
stations  and  along  solid  walls.  Far  field  and  downstream 
conditions  are  remote  from  the  base  and  do  not  affect 
the  solution  in  the  base  region.  Inflow  conditions  for  the 
external  flow  are  supplied  from  experimental  data  at  a 
measurement  station  upstream  of  the  base.  The 
freestream  and  much  of  the  external  boundary  layer  flow 
are  supersonic,  as  is  the  jet  exit  flow.  Therefore,  all 
inflow  variables  are  prescribed.  The  external  boundary 
layer  is  approximated  as  a  l/7th  power  law  mean 
velocity  profile,  and  the  density  and  turbulence  quantities 
are  consistent  with  an  adiabatic,  zero  pressure  gradient 
boundary  layer. 

At  all  solid  walls  the  viscous  no-slip  condition  is 
applied  to  mean  flow  velocities,  the  pressure  is 
extrapolated  linearly,  and  there  is  a  zero  temperature 
gradient.  The  shear  stress  at  the  wall  and  k  and  e  at 
the  first  point  away  from  the  wall  are  determined  by 
incompressible  "wall  function”  boundary  conditions,  which 
are  based  on  the  assumptions  that  the  velocity  profile  is 
logarithmic  and  that  the  production  of  turbulence  kinetic 
energy  equals  its  dissipation.  For  the  flows  considered 


here,  vvall  function  boundary  conditions  are  as  accurate 
as  can  be  obtained  with  any  approach.  Furthermore,  in 
the  base  region,  turbulence  mechanisms  in  the  free  shear 
layers  are  dominant,  and  the  turbulence  model  boundary 
conditions  should  have  a  weak  effect  on  the  solution. 

Grid  Generation 


It  has  become  clear  that  some  form  of  solution 
adaptive  grid  clustering  is  desirable  in  order  to  obtain 
accurate  solutions  affordably.  Deiwert  et  al.'  have 
developed  a  rather  general  method  for  clustering  grid 
points  in  regions  of  high  gradients,  and  they  have  applied 
the  method  to  base  flows.  Venkatapathy  and  Lombard' 
have  used  this  procedure  in  an  approach  that  employs 
adaptive  and  overlapping  grids  for  base  flows.  Despite 
these  efforts,  grid  independent  solutions  are  very  elusive. 

Accurate  solutions  to  the  governing  equations  are 
obtained  by  providing  grid  clustering  which  is  sufficient 
to  eliminate  significant  truncation  errors  at  flow  features 
which  must  be  resolved.  In  the  present  work  it  is 
observed  that  shocks  require  minimal  clustering  because 
their  internal  structures  are  irrelevant  for  the  present 
purposes  and  because  the  differencing  scheme  has  good 
shock  capturing  abilities.  However,  free  shear  layers  (and 
perhaps  expansions)  require  significant  grid  clustering 
because  their  internal  structures  must  be  resolved. 
Denser  grid  clustering  may  be  required  at  free  shear 
layers  for  the  k-e  model  than  for  zero-equation  turbulence 
models  (e.g.,  Baldwin-Lomax),  since  k  and  e  vary  more 
rapidly  through  the  layer  than  do  mean  flow  quantities. 
Therefore,  the  objective  of  the  present  approach  is  to 
provide  the  greatest  clustering  for  the  free  shear  layers 
bounding  the  recirculating  flow,  which  also  gives  good 
clustering  at  the  expansion  corners.  No  special  clustering 
is  done  at  si  'cks.  By  focussing  only  on  the  shear 
layers,  it  is  possible  to  generate  grids  without  excessive 
skewing  or  rapid  variations  in  clustering,  which  can 
degrade  accuracy. 

A  grid  generation  procedure  has  been  developed  that 
makes  use  of  the  fact  that  the  crucial  shear  layers  In  the 
base  flows  considered  here  have  a  relatively  simple 
geometry.  The  shear  layers  emanating  from  the  jet  exit 
and  the  shoulder  of  the  body  are  initially  very  thin  and 
relatively  straight.  They  form  two  sides  of  the  triangular 
shaped  base  region.  (This  approach  is  similar  to 
concepts  in  component  models  for  base  flows.)  The  grid 
generation  method  specifies  "control  surfaces."  which 
follow  the  shear  layers  as  they  depart  from  the  base,  as 
sketched  in  Fig.  3(a).  At  the  control  surfaces  the  radial 
grid  clustering  is  specified,  and  exponential  stretching  is 
used  to  expand  the  clustering  away  from  the  control 
surfaces. 

Where  the  layers  merge,  the  control  surfaces  are 
forced  to  maintain  separation  by  requirements  imposed  on 
the  radial  grid  spacing.  Downstream  of  the  merging  zone 
the  grid  relaxes  to  roughly  uniform  spacing  in  the  radial 
direction,  within  the  plume.  Clustering  is  used  in  the 
axial  direction  near  the  base.  The  axial  grid  spacing  is 
kept  moderately  small  until  somewhat  downstream  of  the 
Mach  disk,  at  which  location  the  flow  is  entirely 
supersonic.  The  minimum  grid  spacing  on  the  base  and 
the  body  is  selected  to  achieve  yn,jn=50,  for  application 
of  wall  function  boundary  conditions. 


3 


fww  pw  V  iwir 


Figure  3(b)  gives  the  grid  in  the  base  region  for  the 
specific  flow  displayed  in  Figs.  4-6,  Looking  ahead  to 
Fig.  5.  it  can  be  seen  that  the  grid  clustering  follows  the 
shear  layers  of  the  solution.  Most  grids  used  in  the 
present  work  are  very  similar  to  that  shown  in  Fig.  3(b). 
For  grid  refinement  studies,  which  are  given  later,  the 
total  number  of  grid  points,  and  thus  the  absolute  grid 
density,  is  changed,  but  the  character  of  the  grid  is 
generally  similar. 


flow  (after  two  grid  refinements).  Of  the  many  solutions 
computed  for  this  flow,  this  one  is  believed  to  be  the 
best  approximation  to  the  solution  of  the  PDFs  which 
model  the  flow  because  it  is  on  the  finest  grid.  The 
present  solution  is  weakly  grid  dependent,  which  will  be 
discussed  later.  The  initial  calculation  on  a  coarser  grid 
gives  a  better  prediction  of  the  base  pressure  (shown  in 
Fig.  13),  but  it  is  believed  to  be  a  less  accurate  solution 
of  the  model  PDFs. 


At  present  the  grid  generation  procedure  is  not 
automated,  nor  is  it  applied  concurrently  with  the  Navier- 
Stokes  method.  However,  there  are  no  fundamental 
research  issues  blocking  these  advances.  The  present 
procedure  is  used  by  first  obtaining  a  converged  flow 
solution,  and  then  using  that  solution  to  guide  the 
generation  of  a  new  grid.  The  approach  is  inefficient  for 
computational  costs,  but  it  minimizes  development  costs. 
Since  no  research  efforts,  including  the  present  one  or 
Refs.  2  or  7,  have  demonstrated  grid  independent 
sol.itions,  there  seems  little  impetus  to  develop  a 
production-type  p'ocedure  at  this  time. 


Results 


Results  are  presented  to  demonstrate  the  accuracy  of 
the  method  and  to  show  specific  areas  where  more  work 
is  needed.  The  two  major  topics  involve  the  effects  of 
the  grid  and  of  the  turbulence  modeling  on  the  accuracy 
of  the  solutions. 


The  afterbody  experiments  reported  by  Petrie  and 
Walker*  are  used  for  comparison.  These  flows  involve 
an  cylindrical  afterbody  with  a  centered  axisymmetric  Jet 
in  a  M  =  1.4  freestream.  The  external  boundary  layer 
at  a  station  3.2  body  radii  upstream  of  the  base  is  0.14 
body  radii  thick,  with  Reg  =  9.0x10*.  The  propulsive 
jet  emerges  from  a  conical  10*  half-angle  nozzle  which 
has  a  design  exit  Mach  number  of  M  =  2.7:  the  ratio  of 
nozzle  exit  radius  to  body  radius  is  0.2.  Inflow  boundary 
conditions  at  the  nozzle  exit  plane  are  consistent  with 
experimental  stagnation  conditions  and  ideal  inviscid 
conical  flow  in  the  nozzle.  Data  were  reported  at  two 
nozzle  pressure  ratios.  NPR  =  2.15  and  NPR  =  6.44. 
Initially,  some  results  are  given  for  both  cases  to 
demonstrate  the  capabilities  of  the  method.  Then  a 
more  detailed  investigation  of  the  factors  which  influence 
the  accuracy  of  the  solution  are  presented,  for  the 
NPR  =  2.15  case.  Results  for  the  variation  of  the 
spreading  rate  of  the  plane  free  shear  layer  with  respect 
to  Mach  number  are  also  given. 

Figure  4  gives  the  base  pressures  computed  by  the 
present  method  and  other  Navier-Stokes  methods  and 
measured  in  the  experiment,  for  the  NPR  =  2.15  case. 
The  present  result  is  reasonably  accurate,  and  it  predicts 
a  relatively  uniform  base  pressure  which  is  a  significant 
feature  of  the  experimental  results. 


Except  for  the  present  results,  the  results  in  Fig.  4^1 
were  reported  in  Ref.  1.  The  results  in  Fig.  4(bJ®’  ’ 
have  been  presented  since  the  publication  of  Ref.  1.  It 
must  be  noted  that  the  other  calculations  shown  in  Fig. 
4(a)  were  perfor^^  "blind."  That  is,  those  performing 
the  calculations^'  '  were  given  operating  conditions  and 
inflow  boundary  data,  but  other  experimental  data  were 
withheld.  The  blind  test  results  give  a  strong  indication 
that  the  predictive  capability  for  this  flow  is  poor.  The 
present  results  were  obtained  with  prior  knowledge  of 
experimental  results,  but  that  knowledge  was  not  used  to 
improve  the  accuracy  of  the  calculation.  The  present 
solution  was  obtained  on  the  third  calculation  of  this 


A  comparison  between  the  predicted  Mach  number 
contours  and  the  experimental  c^adowgraph  in  Fig.  5 
indicates  that  the  general  features  of  the  flow  are 
captured  accurately.  A  clear  distinction  between  the 
barrel  shock  and  the  jet  shear  layer  is  visible  in  the 
Mach  number  contours  but  not  in  the  experimental 
shadowgraph  or  velocity  field  (shown  below).  The 
location  of  the  Mach  disk  is  well  predicted.  The 
computed  maximum  Mach  number  in  the  jet  is  between 
5.0  and  5.5,  which  is  greater  than  the  measured 
maximum  of  about  M  =  3.5.  It  is  not  known  if  this 
discrepancy  is  significant  since  uncertainty  estimates  for 
the  Mach  number  measurements  are  not  given.  The 
Mach  number  may  be  sensitive  to  changes  in  velocity  in 
this  range,  depending  on  how  the  acoustic  speed  is 
determined. 

Velocity  vectors,  given  in  Fig.  6,  reveal  discrepancies 
between  experiment  and  computation  in  the  recirculating 
flow,  which  are  believed  to  be  significant.  In  the 
computation,  the  reverse  flow  velocity  is  too  large,  the 
outer  shear  layer  spreads  too  rapidly,  and  the  vortical 
nature  of  the  recirculation  is  different.  These 
discrepancies  are  small  for  the  present  case,  but  they 
indicate  the  potential  for  larger  errors  in  predictions  of 
other  cases.  The  recirculating  base  flow  is  driven  by  the 
turbulent  stresses,  and  thus,  these  errors  suggest 
shortcomings  in  the  turbulence  model  or  inadequate 
resolution  in  the  shear  layers. 

Results  for  the  case  with  NPR  =  6.44  are  discussed 
briefly  before  further  comments  on  the  NPR  =  2.15  case 
are  given.  The  base  pressures  from  experiment  and 
previous  blind  predictions*  and  from  the  present  effort  are 
given  in  Fig.  7.  The  experimental  shadowgraph  and  the 
predicted  Mach  number  contours  are  compared  in  Fig.  8. 
The  present  results  give  a  reasonable  prediction  of  the 
average  base  pressure  and  the  overall  structure  of  the 
flow  field. 

However,  the  strong  expansion  at  the  lip  of  the 
nozzle  causes  large  oscillations  in  the  pressure.  A  repeat 
calculation  on  a  different  grid  gave  a  significantly 
different  solution.  The  problem  lies  in  the  method's 
inability  to  capture  the  expansion  at  the  nozzle  exit 
which  is  too  fine  to  be  resolved  on  the  grid  employed. 
The  computational  nozzle  exit  flow  does  not  have  a 
physically  correct  boundary  layer,  but  rather  it  has  a  two 
point  transition  zone  between  the  uniform  conical  jet  flow 
and  the  nozzle  lip.  The  validity  of  this  approximation  is 
not  known.  The  physical  nozzle  boundary  layer  is  very 

thin,  so  that  resolving  it  can  be  a  severe  computational 
burden.  Time  constraints  did  not  permit  a  study  of  the 
problems  at  the  strong  expansion;  therefore,  the 
NPR  =  6.44  case  was  not  investigated  further.  Future 
work  will  address  this  problem.  The  expansion  problem 
is  also  present  in  the  NPR  =  2.15  case,  but  it  is 
confined  to  a  few  grid  points  near  the  nozzle  lip  and  the 
global  solution  is  independent  of  changes  to  the  grid  near 
the  lip.  For  the  NPR  =  2.15  case  the  expansion 
problem  can  be  ignored  while  other  aspects  of  the 
calculation  are  investigated. 
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Effects  of  Grid  Resolution 

A  major  concern  is  to  determine  the  accuracy  with 
which  the  numerical  solution  approximates  the  exact 
solution  of  the  model  POEs.  For  purposes  of  research, 
it  would  be  desirable  to  produce  solutions  to  the  model 
PDEs  which  are  independent  of  the  grid  and  numerical 
discretixation  scheme,  except  perhaps  in  the 
neighborhoods  of  discontinuities.  A  grid  dependence 
study  was  performed  in  an  effort  to  determine  if  the 
present  results  were  primarily  due  to  fortuitous  grid 
choice  or  if  they  are  representative  of  the  grid- 
independent  solution  of  the  model  PDEs.  Fig.  9  gives 
the  base  pressures  obtained  for  the  three  different  grids 
used  in  the  grid  refinement  study  and  an  additional  grid 
which  was  specifically  generated  to  give  poor  clustering  in 
critical  regions.  The  changes  in  base  pressure  are  less 
than  +!-  S%  relative  to  freestream  static  pressure,  but 
grid  independerKe  is  not  achieved.  Furthermore,  the  base 
pressure  does  not  display  a  monotonic  trend  as  the  grid 
is  refined.  Thus,  it  is  not  possible  to  estimate  the  base 
pressure  for  an  asymptotically  fine  grid,  but  also,  there  is 
no  indication  that  the  base  pressure  will  change 
significantly  as  the  grid  is  further  refined. 

A  disconcerting  aspect  of  the  present  results  is  that 
they  disagree  with  a  previous  result’  (contributed  by  the 
Lockheed  group*)  obtained  with  a  k-c  turbulence  model. 
Initial  speculation  was  that  inadequate  grid  clustering  and 
first  order  upwind  differencing  on  k  and  e  used  in  Ref.  4 
were  the  cause.  A  central  difference  calculation  was  run 
on  a  grid  similar  to  that  used  by  the  Lockheed  group 
with  a  large  amount  of  artificial  dissipation  added  to  the 
k  and  c  equations  to  mimic  the  first  order  differencing. 
Strong  diffusive  effects  on  k  and  e  were  observed,  but 
the  resulting  base  pressure  was  not  signiflcantly  closer  to 
the  Lockheed  results.  The  blind  cases  were  run  at 
Mm  =  1.343,  (due  to  a  measurement  error  which  was 
later  corrected),  but  a  calculation  at  Mm  =  1.343  also 
failed  to  produce  better  agreement  with  the  Lockheed 
results.  Hence,  there  is  a  significant  difference  between 
the  present  results  and  a  previous  calculation,  obtained 
with  a  similar  method,  which  cannot  readily  be  explained. 

Additional  calculations  were  run  for  the  NPR  =  2.15 
case,  in  which  central  differencing  with  varying  amounts 
of  artificial  dissipation  was  used.  Central  differencing 
degrades  shock  capturing  but  should  not  affect  the 
accuracy  with  which  smooth  fields,  such  as  the  shear 
layers,  are  computed.  The  base  pressures  for  these  and 
all  other  solutions  mentioned  are  given  in  Fig.  10.  This 
gives  an  estimate  of  the  uncertainty  in  the  base  pressure, 
roughly  +/-  5%.  which  can  be  attributed  to  effects  of 
the  grid  and  the  differencing  scheme.  Comparison  of 
Figs.  4  and  10  suggests  that  effects  of  the  grid  and  the 
inviscid  differerKing  scheme  canrrot  account  for.  differences 
in  base  pressure  between  present  and  other’’  ‘  ’*  results. 

Effects  of  the  Turbulence  Model 

In  the  remainder  of  the  results  section,  the  effect  of 
Mach  number  on  the  turbulence  model  is  considered. 
The  modification  to  the  k-e  model  to  account  for  Mach 
number  effects.  Eq.  (3),  was  Incorporated  in  the  method, 
and  calculations  of  a  planar  free  shear  layer  at  several 
Mach  numbers  and  of  the  NPR  =  2.15  base  flow  were 
performed. 

The  spreading  rate  of  the  planar  free  shear  layer  was 
selected  as  a  test  case  et.Lbe  Stanford  Conference  on 
Complex  Turbulent  Flows.’  There  were  two  sets  of 
results  submitted  to  the  conference.  The  results  from 
the  algebraic-stress  two-equation  k-w  model  used  by 


Wilcox  gave  poor  results  and  showed  no  change  of 
spreading  rate  with  Mach  number.  Bonnet's  RSTM 
showed  good  agreement  with  data  up  to  Mm  -  7. 

The  spreading  rate,  as  defined  in  Ref.  13,  given  by 
the  modified  k-e  model  was  evaluated  by  performing 
calculations  of  a  spatially  evolving  shear  layer  for  a  range 
of  Mach  numbers.  The  inflow  conditions  were  rJ  a  flat 
plate  turbulent  boundary  layer  at  Re^  =  9x10  .  The 
fluid  on  the  low  speed  side  of  the  shear  layer  was  at 
rest,  and  it  had  the  same  pressure  and  stagnation 
temperature  as  the  freestream.  The  spreading  rate  was 
evaluated  at  a  downstream  distance  of  x/dg  =  300  (0g 
is  the  initial  momentum  thickness),  at  which  point  the 
spreading  rate  was  nearly  constant  with  respect  to 
distance.  The  spreading  rate  given  by  the  modified  k-6 
model  is  shown  in  Fig.  11.  The  spreading  rate  at 
M  =  0  is  dff/dx  =  0.10,  which  is  slightly  below  the 
experimental  value  0.115.  (This  was  actually  run  at 
M  =  0.5  with  the  Mach  number  modification  removed.) 
This  underprediction  is  consistent  with  k-c  results 
presented  at  the  Stanford  conference  and  reflects  the  fact 
that  the  k-c  model  coefficients  are  optimized  for  a 
variety  of  flows,  not  just  the  plane  free  shear  layer.  The 
computed  spreading  rate  follows  the  trend  of  the 
experimental  data  and  is  in  good  agreement  at  M  =  5. 
The  model  underestimates  Mach  number  effects  slightly, 
but  in  general,  it  gives  a  good  approximation  to  the 
spreading  rate  of  the  free  shear  layer  for  M  ^  5. 

The  modified  model  was  applied  to  the  prediction  of 
the  NPR  =:  2.15  case,  and  the  base  pressure  is  given  in 
Fig.  12.  With  the  modified  model,  the  shear  stress  and 
the  entrainment  in  the  shear  layers  decrease.  This 
causes  the  predicted  base  pressure  to  rise  by  8%  relative 
to  freestream,  so  that  it  is  in  rather  poor  agreement  with 
experimental  base  pressure.  Other  data  (not  shown) 
indicate  that  the  reverse  flow  region  is  now  too  long, 
and  that  a  subsonic  region  between  the  jet  plume  and 
the  freestream  extends  far  downstream,  which  is  also 
contrary  to  experimental  data.  There  are  several  points 
which  should  ^  learned  from  this  result. 

Comparison  of  Figs.  10  and  12  reveals  that  the 
change  of  base  pressure  caused  by  one  modification  to 
the  turbulence  model  (-1-8%)  is  comparable  to  the  change 
in  pressure  which  was  induced  bv  changes  to  the  grid 
and  the  difference  scheme  {+/-  5%).  The  need  for  other 
modifications  to  the  model  is  justifiable.  Thus,  the 
uncertainty  in  the  present  base  flow  calculations  due  to 
turbulence  modeling  is  at  least  equivalent  to  that  due  to 
numerical  resolution.  It  is  difficult  to  argue  that  the 
Mach  number  modification  is  inappropriate  for  this  flow 
or  that  it  is  grossly  in  error  for  the  plane  shear  layer. 
Thus,  there  must  be  additional  errors  in  the  standard  k-e 
model  that  cancel  the  errors  due  to  Mach  number  effects. 
It  is  unlikely  that  this  cancellation  will  occur  in  a  wide 
range  of  flows.  Thus,  good  results  obtained  with  the 
standard  k-c  model,  or  any  model  which  fails  to  predict 
Mach  number  effects,  should  be  viewed  as  fortuitous. 

As  shown  in  Figs.  4  and  7.  most  prediction  methods 
which  use  the  Baldwin-Lomax  model  underpredict  the 
base  pressure.  If  a  Mach  number  modification  similar  to 
Eq.  (3)  were  applied  to  the  Baldwin-Lomax  model,  then 
the  predicted  base  pressures  would  likely  rise  and.  on 
average,  give  better  agreement  with  experimental  data. 


Summary 

A  Navier-Stokes  prediction  method  which  is 
distinguished  by  an  upwind  characteristic  differencing 
scheme  for  inviscid  fluxes  and  a  k-e  turbulence  model 
was  used  to  investigate  the  accuracy  of  base  flow 
calculations.  Grid-independent  solutions  were  not 
obtained.  However,  grid  dependence  was  shown  to  be 
small  relative  to  that  seen  in  previous  work  and 
comparable  to  the  uncertainty  in  the  solution  due  to  the 
turbulence  model.  Despite  its  many  shortcomings  for 
supersonic  base  flows,  the  standard  high-Reynolds  number 
k-e  model  gave  good  predictions  of  base  pressure  and 
many  oth«  aspects  of  the  flowfield. 

However,  when  the  model  was  modified  to  account 
for  high  Mach  number  effects,  the  accuracy  of  the  base 
flow  calculation  was  degraded.  This  suggests  that 
hidden  error  cancellation  is  present  in  calculations  with 
the  standard  k-e  model.  If  this  is  true,  then  similar 
conclusions  apply  to  other  methods.  Good  results  for 
high  speed  flows  from  any  method  which  does  not 
include  Mach  number  effects  on  turbulence  may  be 
fortuitous. 
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Figures 


Fig.  1.  Sketch  of  base  flow  with  propulsive  jet, 
showing  flow  characteristics  and  nomenclature. 
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Fig.  7.  Comparison  of  predicted  ^base  pressure  with 
experiment  and  previous  calculations,'  NPR  =  6.44. 


Fig.  6.  Comparison  of  computed  and  experimental 
velocity  vectors.  NPR  =  2.15. 
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Fig.  9.  Range  of  base  pressures  predicted  by  refining 
grid  clustering,  for  second  order  upwind  differencing 
scheme,  NPR  =  2.15 


Fig.  11.  Spreading  rate  of  a  plane  free  shear  layer  as  a 
function  of  Mach  number. 
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Fig.  10.  Range  of  base  pressures  predicted  by  varying 
grid  clustering,  differencing  scheme,  and  dissipation 
parameter  used  with  central  differencing.  NPR  =  2.15 
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Computation  on  60  x  76  grid, 
Standard  k-c  model 
Mach  number  modified  k-c  model 
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Fig.  12.  Effects  on  computed  base  pressure  of  Mach 
number  modiFication  to  k-c  model,  coarse  grid  calculation. 
NPR  =  2.15. 
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